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Fast algorithms for computing defects and 
their derivatives in the Regge calculus. 



Any practical attempt to solve the Regge equations, these being a 
large system of non-linear algebraic equations, will almost certainly 
employ a Newton-Raphson like scheme. In such cases it is essential 
that efficient algorithms be used when computing the defect angles 
and their derivatives with respect to the leg- lengths. The purpose of 
this paper is to present details of such an algorithm. 



In its pure form (there are variations) the Regge calculus is a theory of gravity 
in which the spacetime is built from a (possibly infinite) collection of non- 
overlapping flat 4-simplices (plus further low level constraints that ensure 
that the manifold remains 4-dimensional everywhere). This construction 
provides a clear distinction between the topological and metric properties of 
such spacetimes. The topology is encoded in the way the 4-simplices are tied 
together while the metric is expressed as an assignment of lengths to each 
leg of the simplicial lattice. 

The Regge action on a simplicial lattice is defined by 
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Abstract 



1 The Regge calculus 




(i.i) 



where 9 U2 is the defect on a typical triangle 02, A a2 is the area of the triangle 
and the sum includes every triangle in the simplicial lattice. This action is 
a function of the leg lengths L 2 (o~i) and its extremization with respect to a 
typical leg length leads to 



where the sum includes each triangle 02 that contains the leg u\. This set of 
equations, one per leg, are known as the Regge field equations for the lattice. 
Constructing a Regge simplicial lattice will require frequent computations of 
the defects as well as frequent solutions of the above Regge equations, a large 
system of non-linear algebraic equations. These are both nontrivial tasks. 

Most articles on the Regge Calculus are concerned with formal issues, such 
as the nature of its convergence to the continuum [1, 2, 3, 4, 5], its use as a 
tool in quantum gravity [6, 7], the relationship of the Regge equations to the 
ADM equations [8, 9] and so on (for further reading see [10, 11]). 

When people turn to numerical computations in the Regge calculus they are 
usually concerned with matters of existence (that solutions exist), accuracy 
(that the solutions compare well with the continuum spacetime) and con- 
vergence (that a sequence of solutions converge to the expected continuum 
spacetime). These questions are usually posed on sufficiently simple lattices 
that the issue of computational efficiency is not a major concern. But in the 
end, if the Regge calculus is to be of any use in numerical relativity, it must 
be at least competitive with contemporary methods (e.g., finite difference 
and spectral codes) in cases such as the merger of binary black holes. Thus 
issues of computational efficiency will be of paramount importance. It will 
not be sufficient to show that the Regge equations can be solved. What will 
be more important is that the solutions can be found with minimal compu- 
tational effort to allow practical long-time evolutions to be undertaken. 

To the best of the author's knowledge, there are no papers that pose the 
specific question What is the best way to compute the defects?. By best we 
mean minimal computational effort, that the cpu time required to compute 
each defect be as short as possible. The main purpose of this article is to 
address that question. We do not claim that our algorithms are optimal but 
rather that they are a significant improvement on contemporary methods. In 
particular we will present details of an algorithm for computing the deriva- 
tives of the defects at virtually no extra cost above that required to compute 
the defects. 
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This paper is organised as follows. In section (2) we introduce various ele- 
ments such as a coordinate frame (the standard frame) and various vectors 
(e.g., normal vectors) that will be used later in section (4) to compute the 
defects. In section (3) we provide a careful discussion regarding angles in a 
Lorentzian spacetime. Explicit equations for the defects are given in section 
(4), and their derivatives in section (5). Finally, in section (6), we present 
some simple timings for our algorithms against traditional methods. 

Throughout this paper we will refer to simplices in two distinct ways. On 
occasions it will be useful to identify the vertices that comprise a 4-simplex. 
In which case we will write (ijklm) for the 4-simplex built from the vertices 
(i), (j), (k), (I) and (m). On other occasions the vertices will be of lesser 
importance and so we will write 04 for a typical 4— simplex. Similar notation 
will be employed for other n— simplices. 

We will adopt the following rules for coordinates indices, tetrad indices and 
vertex labels. We will use Greek letters exclusively for coordinate indices 
(e.g., for expressions valid only in a specific coordinate frame, in particular 
the standard frame as defined in the following section). However Roman 
letters will be used in two distinct ways and thus we choose to divide the 
Roman alphabet into two groups. The first group, a,b,c - ■ -h, will be used as 
tetrad indices while the second group, i,j,k--- will be used as vertex labels 
(i.e., to distinguish one vertex from another). 

2 The standard frame 

It is rather easy to be seduced by the geometric elegance of the Regge calculus 
into thinking that to maintain the purity of the formulation one must not 
only express but also perform all computations in purely geometric form. As 
always, romance gives way to reality and most people find it far easier to 
adopt local coordinates to simplify the computations. The final results are 
always expressed solely in terms of the scalar data (leg-lengths, angles etc.) 
and thus are independent of the chosen coordinate frame. 

The frame which we are about to introduce, which we will refer to as the 
standard frame, is a trivial generalisation (from 3 to 4 dimensions) of a frame 
commonly used in finite element computations on unstructured tetrahedral 
meshes. This frame has been extensively used by others in the field (see for 
example [8, 12]). 

Consider a typical 4-simplex such as (01234). For the standard frame, set 
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the origin of the coordinates at the vertex (0), align the coordinate axes with 
the four edges and set e a := (0a),a = 1,2,3,4 as coordinate basis vectors. 
Now consider any point x within this 4-simplex. Since the metric inside the 
4-simplex is chosen to be flat we can uniquely express the vector (Ox) as a 
linear combination of the e a , namely 

(Ox) = x% (2.1) 

We take the x M to be the standard coordinates of the point x. Note that these 
coordinates are similar but distinct from another popular choice known as 
barycentric coordinates (see [9, 13, 14]). 

Notice if the point x lies on or inside the 4-simplex then the standard coor- 
dinates x M are subject to the constraints 

< x^ < 1 and < x 1 + x 2 + x 3 + x 4 < 1 

If the point x happens to lie on one of the five faces of (01234) then its 
standard coordinates x^ will be subject to the constraints listed in Table (1) 



Index Face Constraints 



1 


(0123) 


x 4 = 


2 


(0124) 


x 3 = 


3 


(0134) 


x 2 = 


4 


(0234) 


x 1 = 


5 


(1234) 


x 1 + x 2 + x 3 + x 4 = 1 



Table 1: The equations that describe the five faces of the 4-simplex (01234) in 
the standard frame. 



The corresponding constraints for the lower order simplices (e.g., (023)) are 
given by suitable combinations of the entries in that table (e.g., for points on 
the face (023) we take the constraints for (0123) and (0234), i.e., x 1 = and 
x 4 = 0). Finally, we note that the coordinates of the five vertices are simply 
given by 

J " 1 = (2 2) 

1 ;i ' i = 1,2,3,4 v ' 1 



4 



The flat metric for this 4-simplex can be constructed from the leg-lengths by 
choosing a constant symmetric 4x4 matrix g^ u such that 

L\ .'/,„- A./'^ A./'^ i,j = 0, 1, 2, 3, 4 

where Aa;^- := x^ 1 — xtf . This leads to a 10 x 10 system of equations for the 
g^u and the solution is easily seen to be 

9ij = \ {Ll + Ll 3 -Ll) i,j = 1,2,3,4 (2.3) 
provided we take Ly = when i = j. 

This completes the definition of the standard frame for the chosen 4-simplex. 
A similar construction can be applied to any other 4-simplex and it should 
be clear that the standard frames of any pair of adjacent 4-simplices will be 
related by a linear transformation. We will make use of this fact later in 
section (4) when we discuss the computation of defect angles. 



2.1 Volumes 

The 4-volume V4 of a typical 4-simplex is given by 

Vi= / / / / VWl dx 1 dx 2 dx 3 dx 4 



0<x 1 ,x 2 ,x 3 ,x 4 <l 
0<x 1 +x 2 +x 3 +x 4 <l 




g\ I I I I dx 1 dx 2 dx 3 dx 4 



(Xx 1 ,x 2 ,x 3 ,x 4 <l 
0<x 1 +x 2 +x s +x 4 <l 
1 r— 

where g = det(^ v ) in the standard frame. Similar results can be derived for 
any n— simplex. If we denote the n— volume (or measure) of the n— simplex 
by V n then 

Vn = - } V¥\ (2-4) 

TV. 

where g = det^g^) in the standard frame for the n— simplex. 
Another popular expression for the volume of a n— simplex is 

o-n/2 
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where W n is the (n + 2) x (n + 2) Cayley-Menger determinant 
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It is a trivial exercise in linear algebra to show the equivalence of this pair 
of equations for V n . 

There is also a useful relationship between V4 and V3 for each of the five faces 
of the 4-simplex. Let (01234) be a typical 4-simplex and choose one of its 
five faces, say (0123) as a base. Let h be the height of (01234) above the 
base (0123). Then a simple extension of base times height theorem to higher 
dimensions gives 

The height h > can be computed by taking the scalar projection of the edge 
vector (in this case (04)) with the inward pointing unit normal to the face 
(0123). If we denote the edge and normal vectors by e a and n a respectively 
then we obtain 

V 4 = ^\e a n a \V 3 (2.5) 



2.2 Normal vectors 

Each 4-simplex is covered by 5 faces, each of which is a 3-simplex. How do 
we construct the outward pointing unit normal to each face? Recall that if 
a plane is described by the equation n^x^ = constant then its normal vector 
is parallel to n M . Thus, by inspection of table (1) we find that components 
of the five normal vectors are as given in table (2). 

Our next task is to compute each of the rij so that each is an outward 
pointing unit vector. 

If is a unit vector than we must have e(rij) = g^n^n" where e(nj) = ±1 
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Index 


Face 


Normal 


1 


(0123) 


nif, = ni5^ 


2 


(0124) 


n 2 ^ = n 2 5l 


3 


(0134) 


ri3ix = n z 5l 


4 


(0234) 


n 4fl = 


5 


(1234) 


n 5M = n 5 (<SJ + 51 + 51 + 6$ 



Table 2: Normal vectors to the five faces of (01234). The rii are scale factors that 
ensure that each vector is an outward pointing unit-vector. 

according to the signature of nf. This leads to 

e (m) = nlg u e(n 2 ) = njg 33 

4 

= n\g 22 e(n 4 ) = rafo 11 e(n 5 ) = ^ ^ 

Since n 2 > we see that that each e(rij) can be computed according to 
e(ni) = sign (# 44 ) e(n 2 ) = sign (# 33 ) 
e(n 3 ) = sign (g 22 ) e(n 4 ) = sign (g 11 ) e(n 5 ) = sign (5) 
where 5 = v=1 g Mt/ and the function sign (x) is defined by 

{ + 1 for x > 
— 1 for x < 
otherwise 

All that remains is to choose an orientation for each nf . 

Consider for the moment the two sides of the face (0123). Notice that the 
leg (04) and the 4-simplex (01234) lie on the same side of this face. Thus if 
we want to be outward pointing we simply demand that 

e(m)n lM (04)" < 

which leads to 

e{n-i)n\ < 
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Since we already know e{n\) this last inequality allows us to correctly choose 
the sign of n 1 . This same logic can be applied to all but the fifth face (which 
we will deal with shortly). This leads to 

< = -sign(/ 4 )^= (2.6) 



V\9 



^ = -sign(s 33 )^= (2.7) 



VW 22] 





n» = -sign(g 22 )^== (2.8) 



< = -sign(^)4== (2.9) 



For the fifth face (1234) the computations are slightly different. For this face 
the vector (04) is outward pointing, thus to ensure that is also outward 
pointing we require 

e(n 5 )n 5M (04)^ > 

which leads to 

e(n 5 )n 5 > 

and thus 

n<t = (g 1 " + g 2 » + + g*) (2.10) 



Now we shall set about constructing two vectors and as partners to 
the vectors n± and n^. For m![ we will require that it be a unit vector, that 
it be orthogonal to both (012) and and finally that it be oriented to point 
away from (012). Similar conditions will be imposed on m^. 

It is not hard to see that, apart from normalisation and orientation, the 
will be of the form 

mi, = mi - Sf^j (2.11) 

m 2 , = m 2 - 5^ (2.12) 

for some choice of numbers m x and m 2 . From here on in the calculations are 
much like those used for nf . We first impose the normalisation conditions 
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leading to 



/ o 34 o 34 \ 
e(mi) = sign U 33 - 



e{m 2 ) = sign I g 

while for the orientation we require 

e(mi)mi M (03) M > and e{m 2 )m 2fl {04Y > 
which leads to 

-1/2 

mi = e(mi) 
m 2 = e(m 2 ) 
2.3 Dot products 



(2.13) 
(2.14) 



33 ff 3 V 4 
» ^44 

34 34 

44 y y 

9 — is" 



-1/2 



(2.15) 
(2.16) 



From the above definitions it is not hard to compute the following dot- 
products (valid only in the standard frame). 



m^m^ = sign g 



34 34 

,33 y y 
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Id 



m^m 2 ^ 



, 44 y 3 V 4 

sign y - 



y 



33 



sign (y 44 ) 
n^n 2M = sign (# 33 ) 



,34 



m^m 2fl 



sign (sV 4 ) " 
<rn 2jU = - sign ( y 44 ) 



34 



mtn 2 ^ 



- sign (y- 33 ) 



1 - 



33^44 1/2 


y u y 34 


1/2 


g 33gU 




y 34 y 34 


1/2 


^44 





These will be used when we develop explicit equations for the defect angles. 
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3 Angles and boosts 



One of the recurrent tasks in the Regge calculus is that of computing the 
angle between pairs of vectors. This would be relatively straightforward if the 
metric was Euclidean. However, for General Relativity, we require the metric 
to be Lorentzian and this introduces some subtleties in the computations. 

3.1 The Euclidean plane 

The following discussion may seem mundane and trivial but it does set the 
scene for the slightly tricky aspects that we are about to encounter in the 
context of a Lorentzian plane. 

Suppose that at a point in a Euclidean plane we have two unit-vectors. How 
do we measure the angle between these vectors? The simple answer given in 
most elementary textbooks is that it equals the length of the arc of the unit 
circle (centred on the given point) bounded by the tips of the vectors. There 
are of course two arcs that join the tips and thus two possible choices for 
the angle. A specific choice is usually made by other considerations (e.g., we 
could choose the arc that follows a counter clockwise direction, which in turn 
requires a definition of clockwise but for the present audience we will take 
that as understood). 

In the Lorentzian plane the set of unit-vectors at a point does not describe 
a closed path but instead four branches of the unit-hyperbola. This posses 
a problem in extending the above definition to the Lorentzian case, namely, 
that it is not clear how to measure the arc-length between points that lie on 
different branches. Thus to avoid this problem we shall re- visit the above 
definition of angles. 

The set of all unit- vectors at a point in the Euclidean plane describes a unit- 
circle. Thus we can establish a one-to-one mapping between the points on 
the unit-circle and the unit-vectors at a point. A simple way to make this 
mapping explicit is to introduce rectangular Cartesian coordinates. Let 
be a parameter that describes a typical point on the unit circle and let u((f>) 
be the corresponding unit- vector with components (u x ((f)) , u y ((f))) . Then we 
have 

u((f>) = e x ii x (<f)) + e y u y ((f)) = e x cos + e y sin (3-1) 

where < < 2tt. Notice that as the unit-circle is traced out in a counter- 
clockwise direction, increases from to In. The mapping is one-to-one so 
a unique inverse exists, which we will denote by <f>(u). 
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Now let <p(u, v) be the counter clockwise angle from u to v. Then we can 
define 0(w, v) by 

(j)(u,v) := <p(v) - (j){u) 

where <fi on the right hand side is obtained by inverting equations (3.1). Note 
that this definition leads to the following properties for angles, namely, 

• (f)(u, —U) = ±7T. 

• <P(v>,v) = -<p(v,u). 

• (f)(u, v) = (fi(u, w) + (p(w, v) for any u, v and w. 

• <f)(u, v) is invariant with respect to rotations of u and v. 

• If u and v are orthogonal then (p(u,v) = tt/2 or 37r/2. 

• 4>{u, v) is bounded. 

The above definitions work well in a 2-dimensional space but for the Regge 
calculus we will be working with 4-dimensions. How then do we gener- 
alise this definition to higher dimensions? The simple answer is to use dot- 
products. Thus, given a flat Euclidean metric g^ v and a pair of unit-vectors 
u and v, the angle v ) can be computed from 

cos 4>{u, v) = g^uU^v" 

Is this definition equivalent to that given above, equation (3.1), i.e., will 
the two equations yield the same angle <p(u, v) for any pair of unit vectors? 
Clearly the answer is yes and the proof is trivial. Expand both vectors u 
and v on an orthonormal basis. We are free to rotate the basis so that the 
vectors are of the form u = e x and v = e x cos a + e y sin a. Upon substituting 
these into the right hand sides of both definitions we see that their left hand 
sides agree, i.e., they give the same angle, in this case 4>(u, v) = a. 



3.2 The Lorentzian plane 

The preceding discussion is all rather simple and could easily be skipped 
over. However it does serve to motivate how we propose to compute angles 
in a Lorentzian space. 
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Unlike the Euclidean case, the set of all possible unit-vectors at a point 
defines four branches of a unit-hyperbola, which we have drawn in figure (1) 
using the parametric form 



e t sinh + e x cosh in I 

e t cosh + e x sinh in II 

et sinh — e x cosh in III 

et cosh — e x sinh in IV 



u(4>) = e t u\4>) + e x u x (<p) = < 



(3.2) 



V 



where u x and u l are the components of a typical radial unit- vector, w(0) = 
e t M*(0) + e x u x ((p). 

Notice that in proceeding along the four branches in a counter-clockwise 
fashion the angle does not increase monotonically. Indeed in regions II 
and IV the angle decreases from +oo to — oo. Why then have we made this 
somewhat peculiar choice? Simply, it is the only choice (apart from changing 
the sign of on all four branches) that ensures that the angle between a pair 
of vectors is invariant with respect to boost transformations. We will return 
to this point in just a moment. 

We can use the above parametric form to suggest a mapping between points 
on the unit-hyperbola and unit- vectors at a point. In this case though the 
mapping is not one-to-one (to each there are four possible unit- vectors u), 
however, given a vector u there exists exactly one value of and thus it is 
meaningful to define as a function of u. Thus as before we can define the 
angle between any pair of unit-vectors in the Lorentzian plane by way of 



where the right hand side is obtained by inverting equations (3.2). 

We now return to the point raised above concerning the non-monotonic be- 
haviour of 0. Had we chosen a different orientation along the second branch, 
such as 



then the pair of vectors defined by u = et cosh l+e x sinh 1 in region II and v = 
et sinh l + e x cosh 1 in region I would have 0(w) = —1 and 0(f) = +1. Thus u 
and v would be separated by an angle of 2. But we recognise that u and v are 
the result of a boost applied to the pair of vectors e t and e x respectively. And 
as the angle between e t and e x is zero we see that the modification suggested 
above violates the boost invariance of the angle between pairs of vectors. 



4>(u,v) := <f>(v) - 0(w) 



ii(0) = etii*(0) + e x u x (4>) = e t cosh — e x sinh 



in II 
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If now we turn to dot-products then we find from the above parameterisations 
that 



where the subscripts on u and v indicate the quadrant in which the vector 
is defined. Though this accounts for just 4 of the 16 possible combinations 
it is all that we need to compute the angle between any pair of vectors. For 
example, for the vectors Un and wni we can use <fi(un, wm) = <f>(win, vi) — 
4>(uu, vi). This result follows directly from the previous definitions given for 
(f)(u,v). 

The following properties of <p(u, v) are readily obtained from the above defi- 
nition (compare these with the similar properties listed in section (3.1)) 

• <p(u, —u) = 0. 

• (j>(u,v) = -<p(v,u). 

• (j)(u, v) = <fi(u, w) + <f)(w, v) for any u, v and w. 

• (p(u, v) is invariant with respect to boosts of u and v. 

• If u and v are orthogonal (and non-null) then <p(u, v) = 0. 

• (f)(u, v) is unbounded. 

This last case can only arise when one or both of the vectors is null. We will 
deal with this problem by simply excluding the case of null- vectors. Fortu- 
nately for most work in numerical Regge calculus it is exceedingly unlikely 
that we would encounter a null- vector (when computing defects). This is far 
from ideal but for the moment we will push on regardless (while acknowledg- 
ing the limitations of our description). 

4 Defect angles 

The customary approach to computing defects is to use what one might 
loosely call a sum of angles method. This goes as follows. Suppose we have 



cosh0(u/, vi) 
sinh0(u/j-, vj) 
cosh <fi(uui, vi) 
■ sinh0(u7v, vi) 
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a 4-dimensional Euclidean lattice. Pick any triangle 02 in the lattice. This 
triangle is surrounded by a closed loop of 4-simplices. Pick any one of these, 
call it cr 4 . This 4-simplex will contain two tetrahedral faces 031 and C32 
that are themselves also attached to 02. Now measure the angle subtended 
at 02 by this pair of tetrahedra, call it 0i2- Finally, repeat this for each of 
the remaining 4-simplices in the closed loop. Then the defect on 02 is simply 
defined as 9 = 2tc — Y^ <Pi2- The dihedral angles are computed (almost without 
exception) by way of 

sin 012 = - (4.1) 
<J ^3,1^3,2 

where V^i and are the volumes of the two tetrahedral faces of 04 while 
V4 is the 4-volume of 04 and V2 is the area of 02 (all of which can be readily 
computed using the equations given in section (2.1)). 

The above prescription is immediately applicable to a Euclidean simplicial 
lattice, however, for a Lorentzian simplicial lattice some minor changes must 
be made. In particular, for a spacelike 02 the right hand side of equation 
(4.1) may have values greater than one. This problem was dealt with by 
Wheeler [15] and Sorkin [16] by allowing the dihedral angles to be complex 
valued. This in turn required that the areas of the triangles also be complex 
valued in order that the Regge action would remain real valued. Though this 
approach is correct it does require some care in choosing the correct branch 
for the inverse sine function (again, see Sorkin [16] for full details). 

We shall take a different (but equivalent) approach to Wheeler and Sorkin 
by working entirely with real valued expressions. This is a minor advantage 
in any computer code as we would only need to reserve half the memory 
required in the Wheeler and Sorkin approach (the defects, in contrast to the 
dihedral angles, are always either purely real or purely imaginary and storing 
such numbers as complex numbers is excessive). 

In the following sections we will provide full details of two methods for com- 
puting the defects, the first uses the sum of angles method while the second 
uses a parallel transport method. 

But before doing so we shall take the opportunity to make a small com- 
ment regarding a common pictorial representation of defect angles in the 
Regge calculus. We will assume for this discussion that we have a simple 
2-dimensional lattice consisting of a set of triangles attached to one vertex. 
A popular representation of such a lattice is to make a cut along one of the 
radial edges so that triangles can be laid out flat in a 2-dimensional space. A 
typical example, for a Euclidean lattice, is shown in figure (3). The shaded 
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wedge is the part of the space not covered by the lattice. Also, to recover 
the lattice, we must identify the radial edges of the wedge. In this picture 
we have also included the unit circle so that we can easily measure the angle 
subtended by the wedge and clearly this angle is exactly equal to the de- 
fect angle of the lattice. It should also be clear that the cut can be made 
along any radial line starting from the bone. Choosing some other radial 
leg produces a similar diagram but with the wedge rotated about the bone. 
The flexibility to locate the cut anywhere in the lattice does not apply for a 
Lorentzian lattice. This is easily seen by reference to figures (4-6). In figure 
(4) we see that the defect angle would be positive yet in figure (5) we have 
a negative defect (the angle decreases along this branch when traversed in 
the counter-clockwise direction). Finally, in figure (6) we show an impossible 
construction - the edges of the wedge can not be identified because their 
tangent vectors have differing signatures. 



4.1 Sum of angles 

The definition here is very simple. The defect 9 U2 on a typical bone a 2 is 
defined by 

2n — (0i2) CT4 f° r cr 2 timelike 

0-4(0-2) 

(4.2) 

- 2J (0i2),r 4 f° r a 2 spacelike 

04(02) 

where the sum includes all of the 4-simplices a 4 attached to the bone a 2 . 
The dihedral angles 0X2 can be computed using equation (4.1) taking due 
care for the complex numbers that might arise. On the other hand we can 
work solely with real numbers by combining our definitions of angles given in 
section (3) together with the dot-products of the normal and tangent vectors 
(as defined in section (2.2)). Then for a timelike bone we find 

012 = arccos (m^m 2fl ) (4.3) 

while for a spacelike bone we find 

012 = sign (mlmifj) arcsinh p 12 (4-4) 

where pi 2 is computed from 

{sign (nim 2l /) when \m^m 2tJ \ < |nj t m 2M | 
(4.5) 
sign {m\m 2v ) n^m 2 ^ in all other cases 

Recall that the dot-products in the standard frame are given in section (2.3). 
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4.2 Parallel transport 



Here we will use a commonly quoted but rarely used method for computing 
a defect angle. After parallel transporting a vector around a simple closed 
loop the defect angle is computed as the rotation angle (or boost) between 
the initial and final vectors. If we denote the initial and final vectors by v 
and v' respectively, then we expect a transformation of the form 

v = Rv 

where R is a 4 x 4 rotation (or boost) matrix (this matrix will act trivially 
on vectors parallel to the bone). 

We will construct R by taking account of the separate transformations that 
arise, first from the parallel transport within a 4-simplex, and second from 
the passage from one 4-simplex to the next. In each case we will employ 
suitably chosen basis vectors local to each 4-simplex. The successive stages 
of the transformation will be recorded as a matrix of scalars with respect to 
the local basis. 

We begin by considering the transformation associated with the parallel 
transport within a typical 4-simplex, (ijkl2). The section of the closed loop 
that passes through this 4-simplex will have entry and exit points on distinct 1 
faces of that 4-simplex. Suppose that the entry face is (ijkl) and the exit 
face is (ijk2). At the entry point we will use the (non-orthogonal) tetrad 

ei = (ij) e 2 = (ik) e 3 = n^e^ e 4 = m^e M 

while at the exit we will use 

e'x = (ij) e 2 = (ik) e' 3 = ra£e„ e' 4 = m£e M 

where the e M are the coordinate basis vectors in the standard frame and n^, 
rig, 77ij* and m2 are defined in section (2.2). Then the parallel transport of 
the vector v along the path within this 4-simplex is described by 

v = v a e a = v' a e' a 

where the v a and v' a are the components of v onto the respective tetrads. 
We also know that the entry and exit tetrads are related by a simple linear 
transformation of the form 

ei = e\ e 2 = e' 2 e 3 = T|e 3 + T 3 V 4 e 4 = T 4 3 e' 3 + T 4 4 e 4 (4.6) 

: We exclude the case where the entry and exit faces are one and the same, e.g., when 
the path doubles back on itself. 
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for some set of numbers (which we will compute explicitly in the fol- 
lowing section). This in turn implies a similar linear transformation on the 
components, v a and v a , 



,'3 



T> 3 + TjV 



Ttv 3 + TV 



This completes the transport of the vector within one 4-simplex. Now we 
must take account of the changes that occur as the vector is transported 
across the junction between a pair of 4-simplices. This is very easy to do, we 
need only recognise that the normal vector of the exit face (of the current 
4-simplex) is oppositely oriented to the normal vector on the entry face (of 
the next 4-simplex). This amounts to the following transformation 



v a i-» S$v b 



where S = diag(l, 1, —1, 1). 



We are now ready to compute the transformation matrix R for the closed 
path around the bone. Let T m be the matrix built from the T\ associated 
with the m th 4-simplex, then we have 



R=(ST n ) (5T n _ 1 )(5T n _ 2 )-..(5T 1 ) 



(4.7) 



What can we say about the form of R7 We know that after one loop of 
the bone any vector initially parallel to the bone will be returned unchanged 
while any vector orthogonal to the bone will be subjected to a rotation (for 
a timelike bone) or a boost (for a spacelike bone). Thus we expect R to be 
of the form 



R 



and 



R 



1 














1 














cos (3 


— sin (3 








sin/3 


cos/3 


1 














1 














cosh (3 


sinh (3 








sinh/3 


cosh f3 



for a timelike bone 



(4.8) 



for a spacelike bone 



(4.9) 



for some number (3. The question now is, what is the relationship between (3 
and the defect angle 91 We will answer this question by requiring that the 
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defect angle calculated by the parallel transport method agrees with that 
given by the standard sum of angles method. 

We begin with the simple case of a timelike bone with a defect angle 9. To 
bring the sum of angles method into the picture we will use the representa- 
tions introduced in section (4). Thus in figure (8) we have a timelike bone 
with a defect angle 9. We have also included in that figure the results of 
parallel transporting a pair of vectors n, m along a closed path around the 
bone. Note that the shaded region is excluded from the figure and that the 
radial lines should be identified. We have also suppressed the two dimen- 
sions parallel to the bone (there is no useful information contained in those 
dimensions). Equally, the figure could be viewed as an isometric mapping 
of the 2-dimensional subspace orthogonal to the bone into a subset of R 2 
with a Euclidean metric. After one loop of the bone we obtain the pair of 
vectors n', m' and these are to be compared with the original vectors n", m". 
Evidently this entails a rotation through an angle /3 and so we deduce that 
the defect is given by 9 = (3. 

A little more care is required for the case of a spacelike bone. As we shall 
soon see, the entries in the matrix R may depend upon which 4-simplex is 
chosen to start the loop around the bone. Consider the example shown in 
figure (9). This is very similar to the previous figure with two important 
differences, first we are now working with a Lorentzian metric (in this 2- 
space) and second we must now speak of boost parameters (as opposed to 
rotation angles). Thus the skewed appearance of n, m in the figure is simply a 
consequence of using a Lorentzian metric. By inspection we can see that the 
boost parameter (for the boost from n", m" to n', m') is negative while the 
defect is positive. Thus in this instance we have 9 = —f3. Now turn to figure 
(10). This too is also for the case of a spacelike bone but here we have made 
the cut at a position which has n as a spacelike vector (in the previous case 
n was a timelike vector). Again we see that the boost parameter is negative 
but on this occasion we have a negative defect, thus we have 9 = +/3. 

This result can also be understood by purely algebraic means. Consider the 
transformation matrices that would arise by completing the parallel trans- 
port around a loop by starting in distinct 4-simplices. We expect that the 
defects computed from these matrices must all agree (they describe the paral- 
lel transport around the same bone). In figure (7) we have drawn the first few 
4-simplices attached to the bone (ijk). On each 4-simplex we have drawn 
their particular entry basis vectors (and for simplicity we have not drawn 
them as skewed vectors). Let R\ be the parallel transport matrix associated 
with (ijkl2), with R 2 for (ijk23) and so on. Suppose that n\ and ri2 are 
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timelike while n§ is spacelike. Let P12 be the transformation matrix that 
maps ni,rni into n 2 ,m 2 . Then it follows that 

Rl = P l2 R 2 P\ 2 



However, since n\ and n 2 are both timelike we see that P\ 2 must be a pure 
boost. Furthermore, P i2 and P 2 act in the same 2-space and thus they 
commute, so we have 



Pi — -Pi2 1 -R2-Pi2 — P\2 P12R2 



R? 



which leads to Pi = f3 2 . If we attempt the same calculation for the simplices 
(ijkl2) and (ijk56) we must take account of the fact that n\ is timelike while 
n 5 is spacelike. The map P15 from ni,mi to ns,m 5 will now be of the form 



Pi 



15 



B15Q 



where P15 is a pure boost and (if we ignore the two dimensions parallel to 
the bone) 

" 1" 
-1 



Q 



This leads to 



Pi — 

= Q^ 1 Pfg 1 R 5 B 15 Q 

= Q-'RsQ 



and thus fii 



"As 



So we now have a simple answer to the earlier question, how is 9 related to 







— arcsin (P3 4 ) for a timelike bone 
e(ni) arcsinh (P3 4 ) for a spacelike bone 



(4.10) 



4.2.1 The standard frame 

We can compute the T\ entries in each transformation matrix T m by taking 
appropriate scalar products of (4.6). This leads to 

T 3 = m^2 T 4 = n ltl m% 

3 „ „ 7/ 3 



n 2u n 2 m 2u m 2 
3 _ m^rag . _ 

n 2v n v 2 m 2v m v 2 
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with all other T\ = except = T| = 1. In the standard frame we find 
that 



rp3 
1 3 



T 



T 



T: 



sign (,rV ; </ :l V ; 1 

-sign(^)/ 4 |^VT 1/2 

1/2 



1 



9 M 9 M 



g Z3gU 



sign I 1 



9 U 9 U 

g 33gM 



1 - 



9 34 9 U 



g 33gU 



1/2 



(4.11) 
(4.12) 

(4.13) 
(4.14) 



From the above it is not hard to show that det (ST m ) = 1. 

In summary, we would use equations (4.11-4.14) to compute the successive 
T m , these would then be used to compute R from (4.7) and finally equation 
(4.10) would be used to compute the defect. 

There is one further subtlety in the transformation process that we neglected 
to mention. For a general path (not necessarily the closed path around a 
bone) the 2-simplex that lies at the intersection of the entry and exit faces 
will change along the path. This has not been accounted for in the above 
analysis but the adaptions required are simple and follow along lines similar 
to that given above. Consider two 4-simplices (ijkl2) and (ijk23). Suppose 
the entry and exit faces for the first 4-simplex are (ijkl) and (ijk2) while in 
the second 4-simplex the entry and exit faces are (ijk2) and (jk23). Thus 
upon entry to the second 4-simplex the tetrad will be 

d = (ij) e 2 = (ik) e 3 = e 4 = m^e M 



while at the exit we might use 



e'i = (jk) e' 2 = (j2) e' 3 = <e M e 4 = m^e M 

and this will impose a further transformation on the tetrad components prior 
to leaving the second simplex. We choose not to give the full details here as 
they are not required for the case of most interest, namely, the computation 
of the defect angle. 
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5 Derivatives of defects 



One particularly attractive feature of the standard frame is that it allows us 
to compute the derivatives of the defects at virtually no extra computational 
cost beyond that required for the defects. In short, we get the derivatives for 
free. 

We begin by considering one 4-simplex <7 4 := (ijkl2) and focusing on the an- 
gle 0i2 between the pair of 3-simplices 03(1) := (ijkl) and 03(2) := (ijk2). 
Suppose for the moment that a 2 '■= [ijk) is timelike then with the orienta- 
tions for n M and m M chosen as per figure (2) we have 

m 2fl = -n lfl sin i2 + m lM cos i2 
n 2(l = -n lfJ , cos 0i2 - m Xtt sin <f> 12 

If some small changes are now made to the lattice (e.g., by small changes in 
the leg-lengths) then we must have 

5m 2 ^ = n 2fM S(pi2 - Sn lfl sin 12 + 5m lfl cos <f> 12 
Sn 2f , = -m 2At 50i2 - 5n lfl cos <f) 12 - 5m ljU sin0i 2 

and thus 

250i 2 = ni^dm^ - m 2 bn 2ill + m^bu^ + n^5m llM 

This result applies in all frames but when restricted to the standard frame 
it takes on a particularly simple form. In the standard frame we have 

m 2il = m 2 (8* + a 2 Sl) n 2il = n 2 <^ 

where the m« and rii are normalisation factors while a, is chosen so that 
= rii^m^. Thus we have 

= m^5n iv for i — 1,2 

while from = n^mf , = n^m ifl and = g^n^m^ we also have 

nfSrriif, = Sg^nfm? for i = 1, 2 

which leads to the following simple equation for the variations in the angle 

25012 = (^2 n 2 + m l n l) &9nv 
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An almost identical analysis can be applied in the case of a spacelike bone. 
It begins with 

irij = — sinh 0i 2 + cosh 0i 2 
i%2 = —n^ cosh 12 + sinh 12 

and leads to 

25012 = - («l2 n 2 + 

To compute the derivative of the defect 9 on <7 2 we need only combine the 
above results with the equation (4.2) for the defect to obtain 

= h (CT2) E ( (m ^ + mW) t^) (5 - 1} 

^ ' CT 2 CT4 ( CT2 ) ^ ' CT 4 

where the sum includes each of the 4-simplices attached to 02, L 2 is a typical 
(squared) leg-length and tip-i) = —1 for a timelike bone and +1 for a spacelike 
bone. 

The utility of this equation should not be overlooked. The partial derivatives 
of the metric in the standard frame are simple numbers such as 1, ±1/2 and 
zero and thus there is no extra cost (of any importance) in computing the 
Jacobian of the Regge equations in situ with the equations themselves. All of 
the terms in this equation are already in use during the computation of the 
defects. This is a significant advantage when it comes to solving the Regge 
equations by standard numerical methods, having the Jacobian at hand at 
no cost is a great bonus. Of course the Jacobian could be computed by finite 
differences but that would be extremely expensive. Each leg would need to 
be varied in turn and the corresponding changes in the defects recorded. If 
we suppose that, on average, each defect depends on N legs, then this finite- 
difference process will increase the computational cost of the defects by a 
factor of N, a significant expense (a bare minimum for N is 10, the 10 legs 
of one 4-simplex). 



6 Timing 

The stated aim of this paper is to present efficient algorithms for computing 
the defects and their derivatives. Here we will present the results of some 
simple tests. We chose two simple lattices each consisting of one bone, a 
timelike bone in one lattice and a spacelike bone in the other. The leg 
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lengths were assigned by choosing a well known metric (e.g., Schwarzschild, 
plane waves, Kasner etc.) and using a geodesic integrator to compute the 
lengths of short geodesies (see [4] for full details). The results presented 
here are a measure of the arithmetic complexity of each algorithm. Thus it 
is reasonable to expect that the our results are largely independent of the 
details of the lattice (e.g., the size, number and location of each 4-simplex in 
the parent spacetime). 

In our first test we computed the defects using the sum of the angles method 
as well as the parallel transport method. The angles required in the sum 
of angles method can be computed in two ways, by the computing ratios of 
volumes or by dot products of normal vectors. In table (3) we have listed the 
cpu times to compute the defects for each of these methods. The times are 
expressed as ratios relative to the standard algorithm (i.e., where angles are 
computed by ratios of volumes) . This shows that the dot product and parallel 



Algorithm 


Equations 


CPU time 


Ratios of volumes 


(4.1) 


1.00 


Dot products 


(4.3,4.4) 


0.58 


Parallel transport 


(4.10) 


0.57 



Table 3: CPU times for three algorithms used to compute the defects. 

transport algorithms are comparable in speed and that both are nearly twice 
as fast as the standard algorithm. This is a good start. 

Now we turn to the computation of the derivatives of the defects. Here 
we will work solely with the dot product algorithm (the angles by ratios of 
volumes algorithm is too slow while the parallel transport method runs neck 
and neck with the dot product method). 

In this test we compute the cpu time to compute the derivatives of the 
defects with respect to the three legs of the bone. We do so using two 
algorithms, the first uses the exact expression (5.1) derived in section (5) and 
the second uses centred finite differences on equations (4.3,4.4). The results, 
normalised against the cpu time to compute the defect, are listed in table 
(4). This clearly shows that there is very little overhead when computing the 
derivatives using the exact equations (5.1). In contrast the finite difference 
method is close to six times slower. This figure is easily understood - for 
each of the three legs we have to compute the defect twice, hence a total of 
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Algorithm Equations CPU time 



Exact derivatives 



Finite differences 



(5.1) 
(4.3,4.4) 



6.03 



1.11 



Table 4: CPU times to compute the three derivatives of the defect. 

six defect computations. Now consider a typical bone in a typical simplicial 
lattice. The defect on this bone will be depend on large number of legs, at 
least 10 and possibly upwards of 50. The above calculations suggest that 
computing all the derivatives in this case could be as bad as 100 times slower 
than by the exact equations. This would be so slow as to be impractical 
(assuming no other part of the code dominated the computations). 

7 Cross checking 

It is a common (and regrettable) fact that errors of many kinds (numerical, 
coding, logical) do creep into computer programs. The wise programmer will 
employ as many tests as they can to validate their code. For us we have two 
identities that can be usefully employed, Stokes' theorem for one 4-simplex, 



which can be used to test that the unit normals are correctly oriented and 
Regge's identity 



which can be used to check the correctness of the defect angles (by employing 
2nd order finite differences and looking for 2nd order convergence of the right 
hand side to zero). 



5 




(7.1) 



i=i 




(7.2) 
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7.1 Stokes' theorem 



A typical 4-simplex, say (01234), has five faces. In section (2.2) we proposed 
a simple form for the unit normal to each face from which it follows that 

4 

b i=i 1 

Out of this equation we can now construct a Stokes' theorem for one 4- 
simplex. Consider the first face, say (0123), and its outward normal vector 
in the form ni^ = n^. Upon substituting this into equation (2.5) we find 

4U 4 = \ni\V 31 = e(ni)niV3i 

where V4 is the 4- volume of the 4-simplex and V31 is the 3-volume of the first 
face. Clearly we can repeat this for each of the remaining faces leading to 

4U 4 = e(ni)niV 31 = e(n 2 )n 2 V 3 2 = e(n 3 )n 3 V 33 = £(^4)^4^34 = -e(n 5 )n 5 V 35 

This can be used to eliminate the in the first equation of this section, the 
result is 



= ^ e(ni)n itl V3i 
i=i 



7.2 Regge's identity 

Here we will use equation (5.1) for the derivatives of the defect to obtain a 
simple proof of Regge's identity (7.2). We being by forming a weighted sum 
of equation (5.1) over all of the bones in the interior of the lattice 

<T 2 V / °"2 a-2 0-4(0-2) V / °" 4 

This sum contains contributions from each 4-simplex and its component sub- 
simplices (in particular the faces and the bones). Upon careful inspection we 
see that we can re-order the sums as follows 




Finally, by applying Stokes' theorem to each a 3 we see that the innermost 
sum vanishes and thus we recover Regge's identity (7.2). 
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8 Historical note 



The seeds of the theory that would later be known as the Regge calculus 
were sown in a barber's shop in Princeton sometime in 1959 [17]. Tullio 
Regge, then a PhD student working with John Wheeler, sat in the barber's 
chair and stared into the mirrors that covered the four walls and ceiling of 
the shop. In these mirrors he saw the barber shop repeated over and over 
and this led him to muse, while the barber went about his business, that 
perhaps the universe might be built from, or approximated by, a large array 
of interconnected cells. In 1961 Regge published his seminal paper, General 
Relativity without coordinates. The subject was championed by John Wheeler 
in Groups Relativity and Topology, who coined the title Regge calculus. 
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Figure 1: This figure shows our choice of angle in a 2-dimensional Lorentzian 
space. Note that the angle does not increase monotonically as we traverse the 
branches in a counter-clockwise direction. 
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Figure 2: A typical 4-simplex (ijkl2). Here we show the outward pointing unit- 
vectors ni, ri2 and the corresponding unit tangent vectors mi and m2- These 
vectors are used to compute the dihedral angle 4>i2- They are also used as basis 
vectors for the parallel transport of vectors around the bone (ijk). The small circle 
is our way of noting that (ijk) is a 2-simplex (whereas vertices are drawn as a solid 
dot). 
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Figure 3: An example of a 2-dimensional lattice (with just one bone) cut open 
and laid flat in a Euclidean space. The shaded wedge is not part of the lattice and 
its radial edges should be identified. The defect for this bone is exactly equal to 
the angle subtended by the wedge. The unit circle is drawn simply to remind us 
that the signature is Euclidean. 
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Figure 4: A similar construction to figure (3) but in this instance for a Lorentzian 
signature. Notice the four branches of the unit hyperbola, these match those from 
figure (1). In this example the defect is positive. 
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Figure 5: Here we have a case were the triangles were cut along a timelike edge. 
In this case we have a negative defect angle (and thus this is not the same lattice 
as shown in the previous figure). 
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Figure 6: This arrangement is impossible. The two radial edges of the wedge 
have different signatures and thus they can not be identified. 
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Figure 7: Here we show the first few 4-simplices attached to the bone (ijk). Each 
pair of vectors here are the basis vectors for the entry face for each 4-simplex. This 
figure is not a map of the simplices into a single flat space in the manner used in 
figures (4-6). It is just a convenient picture of the vectors and their relationships 
to their simplices. 



33 




Figure 8: Here we show how the vectors n and m may be parallel transported 
around a bone. In contrast to figure (7) here we do map the 4-simplices into a 
single flat space, hence the introduction of the shaded wedge. This example is for 
a timelike bone. By inspection we see that 9 = (3. 
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Figure 9: This is similar to the previous figure but modified for the case of a 
spacelike bone. We are now working with boosts rather than rotations. Here we 
assume that m is spacelike and thus we must have j3 < 0. We also see that > 
thus we have 6 = —(3. This is the situation we expect if we do a parallel transport 
starting from (ijkl2) (see figure 7) 
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Figure 10: This is also for a spacelike bone but now our parallel transport loop 
starts and finishes in (ijk56). Here we see that < and now < (again, this 
is not the same lattice as used in the previous figure). Now we find 9 = (5. 
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